#include <sstream>
#include <filesystem>
#include <iostream>

#include "gdal/gdal.h"
#include "gdal/gdal_priv.h"

#include "LibGridTIFF/GridTIFF.hpp"
#include "LibDGGSKernel/IDGGSKernel.hpp"
#include "CuLibDGGSKernel/ICuDGGSKernel.cuh"

#include "Common.cuh"
#include "Ephemeris.hpp"

namespace fs = std::filesystem;
using namespace NNU::DGGS;
using namespace NNU::DGGS::Base;
using namespace NNU::DGGS::Kernel;

/// 输入参数
int pGridDimX;
int pBlockDimX;
DGGSGridType pGridType=DGGSGridType::Unknown;
DGGSElementType pElementType=DGGSElementType::Unknown;

///全局变量，减少函数调用开销
ICuDGGSKernel*pCuDGGSKernel=nullptr;

int globalGridSize;                                  //GridTIFF大小
int globalOriginLevel;                               //初始层级
int pLevel= 0;                                       //数据层级
DGGSDataType globalDataType;                         //数据类型

DGGSGridTIFFCode globalGridCode;
DGGSGridTIFFCode *bufferGridTIFFCodes = nullptr;    //缓冲区GridTIFF
int bufferGridTIFFCount;                            //缓冲区GridTIFF数量
void *bufferGridTIFFValues= nullptr;                //缓冲区高程值
double bufferBound[4];                              //缓冲区包围盒


#pragma region device function Declaration

/// <summary>
/// 根据经纬度、坐标方位角、距离求取下一个点的经纬度
/// </summary>
/// <param name="lon">已有点经度</param>
/// <param name="lat">已有点纬度</param>
/// <param name="yaw">坐标方位角</param>
/// <param name="dis">距离</param>
/// <param name="lon1">待求点经度</param>
/// <param name="lat1">待求点经度</param>
__device__ void getNextCoord(double lon, double lat, double yaw, double dis, double& lon1, double& lat1);


/// <summary>
/// 获取格网值
/// </summary>
/// <param name="row">行号</param>
/// <param name="col">列号</param>
/// <param name="srcValueList">GridTIFF值</param>
/// <param name="srcTrans">GridTIFF转换参数</param>
/// <param name="gridSize">格网大小</param>
template<class T>
__device__ T getValueFromGridTIFF(int row,int col,T*srcValueList,GridTransform srcTrans,int gridSize);


/// <summary>
/// 获取格网所属GridTiff的索引号
/// </summary>
/// <param name="basePartition">基菱号</param>
/// <param name="col">列号</param>
/// <param name="row">行号</param>
/// <param name="gridTiffCodes">GridTIFF组</param>
/// <param name="gridTiffCount">GridTIFF数量</param>
/// <param name="originLevel">初始层级</param>
__device__ int getGridTiffIndex(int basePartition,
                                int row,
                                int col,
                                const DGGSGridTIFFCode* gridTiffCodes,
                                int gridTiffCount,
                                int originLevel
                                );


/// <summary>
/// 获取GridTiff的转换参数
/// </summary>
__device__ GridTransform getGridTrans(DGGSGridTIFFCode gridTiffCode,int originLevel);


/// <summary>
/// 判断卫星是否可见
/// </summary>
/// <param name="satSpaceCoord">卫星坐标</param>
/// <param name="basePartition">格元基菱号</param>
/// <param name="row">行号</param>
/// <param name="col">列号</param>
/// <param name="level">层级</param>
/// <param name="eleValues">高程数据</param>
/// <param name="griTIffCodes">相关GridTiff</param>
/// <param name="griTiffCount">相关GridTiff数量</param>
/// <returns></returns>
template<class T>
__device__ int isVisible(SpatialCoord satSpaceCoord,
                         int basePartition,
                         int row,
                         int col,
                         int level,
                         T*  eleValues,
                         DGGSGridTIFFCode*griTiffCodes,
                         int griTiffCount,
                         int originLevel,
                         ICuDGGSKernel *isea);


/// <summary>
/// 计算PDOP
/// </summary>
/// <param name="satCoords">卫星坐标</param>
/// <param name="satVis">可见卫星id数组</param>
/// <param name="satVisNum">可见卫星数</param>
/// <param name="tagSpaceCoord">地面点坐标</param>
/// <returns></returns>
__device__ float calPdopValue(SpatialCoord* satCoords, int* satVis, const int satVisNum, SpatialCoord tagSpaceCoord);


/// <summary>
/// 卫星可见性矩阵组计算
/// </summary>
/// <param name="satvisinfo_d">卫星可见性矩阵组</param>
/// <param name="d_satTime">卫星星历数据</param>
/// <param name="tinfo_d">六边形格网菱形瓦片数据</param>
/// <param name="rowArray_d">行号数组</param>
/// <param name="colArray_d">列号数组</param>
/// <param name="baseP">基菱形块号</param>
/// <param name="baseRow">菱形瓦片行号</param>
/// <param name="baseCol">菱形瓦片列号</param>
/// <param name="triParaments_d">正二十面体参数</param>
/// <param name="vertTable_d">正二十面体参数</param>
/// <param name="itime">卫星时刻</param>
/// <returns></returns>
template<class T>
__global__ void getSatVisMat(SpatialCoord*satSpatialCoords,
                             int  satIndex,
                             T*   eleValues,
                             DGGSGridTIFFCode* griTIffCodes,
                             int  griTiffCount,
                             int  originLevel,
                             DGGSGridTIFFCode tagGridTiffCode,
                             int* visValues,
                             ICuDGGSKernel *isea
);


/// <summary>
/// PDOP矩阵组计算
/// </summary>
/// <param name="satvisinfo_d">卫星可见性矩阵组</param>
/// <param name="d_satTime">卫星星历数据</param>
/// <param name="tinfo_d">六边形格网菱形瓦片数据</param>
/// <param name="rowArray_d">行号数组</param>
/// <param name="colArray_d">列号数组</param>
/// <param name="resPDOP_d">PDOP解算结果存储数组</param>
/// <param name="baseP">基菱形块号</param>
/// <param name="baseRow">菱形瓦片行号</param>
/// <param name="baseCol">菱形瓦片列号</param>
/// <param name="triParaments_d">正二十面体参数</param>
/// <param name="vertTable_d">正二十面体参数</param>
/// <param name="itime">卫星时刻</param>
/// <returns></returns>
template<class T>
__global__ void getPDOPMat(T*eleValues,
                           int* visValues,
                           SpatialCoord*satSpatialCoords,
                           int  satCount,
                           float*tagPDOPValues,
                           GridTransform tagTrans,
                           int  gridSize,
                           ICuDGGSKernel* isea
);


#pragma endregion

#pragma region device function Implementation

__device__ void getNextCoord(double lon, double lat, double yaw, double dis, double& lon1, double& lat1)
{
    double re = double(earthRadius*1000) * cos(double(lat * degToRad));
    lon1 = lon + dis * sin(double(yaw * degToRad)) / double(re * 2 * M_PI) * 360;
    lat1 = lat + dis * cos(double(yaw * degToRad)) / double(re * 2 * M_PI) * 360;
    if(lon1<-180.0)
        lon1=lon+360;
    if(lon1>180.0)
        lon1=lon-360;

    if(lat<-90.0)
        lat=-90.0;
    if(lat>90.0)
        lat=90.0;
}

template<class T>
__device__ T getValueFromGridTIFF(int row,int col,T*srcValueList,GridTransform srcTrans,int gridSize)
{
    int gridX, gridY;
    code2XY(row, col, srcTrans.rowOrigin, srcTrans.colOrigin, gridSize, gridX, gridY);
    int width = gridSize * 2;
    int height = gridSize * 2;
    T value;
    if (gridX >= 0 && gridX < width && gridY >= 0 && gridY < height - 1)
        value=srcValueList[(height - 2 - gridY) * width + gridX];
    else
        value=(T)noData;
    return value;
}

__device__ int getGridTiffIndex(int basePartition,int row,int col, const DGGSGridTIFFCode* gridTiffCodes,int gridTiffCount,int originLevel)
{
    int level = gridTiffCodes->level;

    long long morton = 0;
    long long tiffRow = row / (long long) pow(2,originLevel);
    long long tiffCol = col / (long long) pow(2,originLevel);
    for (int i = 0; i < level - originLevel; ++i)
    {
        morton |= (tiffCol & (1LL << i)) << i | (tiffRow & (1LL << i)) << (i + 1);
    }
    for(int i=0;i<gridTiffCount;i++)
    {
        if(gridTiffCodes[i].basePartition==basePartition&&gridTiffCodes[i].morton==morton&&gridTiffCodes[i].level==level)
            return i;
    }
    return -1;
}

__device__ GridTransform getGridTrans(DGGSGridTIFFCode gridTiffCode,int originLevel)
{
    long long row = 0;
    long long col = 0;
    int level = gridTiffCode.level;
    long long morton = gridTiffCode.morton;

    for (int i = 0; i < level - originLevel; ++i) {
        row |= (morton & (1LL << (2 * i + 1))) >> (i + 1);
        col |= (morton & (1LL << (2 * i))) >> i;
    }

    row = row * (long long)pow(2, originLevel);
    col = col * (long long)pow(2, originLevel);

    GridTransform trans;
    trans.basePartition= gridTiffCode.basePartition;
    trans.rowOrigin=row;
    trans.colOrigin=col;
    trans.level=gridTiffCode.level;
    return trans;
}

template<class T>
__device__ int isVisible(SpatialCoord satSpaceCoord, int basePartition, int row, int col, int level, T* eleValues, DGGSGridTIFFCode*griTiffCodes, int griTiffCount, int originLevel, ICuDGGSKernel *isea)
{
    int gridSize=(int)pow(2,originLevel);

    DGGSElementCode code;
    code.level=level;
    code.elementType=DGGSElementType::Cell;
    GeoCoord coord;

    //格元经纬度坐标
    double originLon=0, originLat=0;
    code.basePartition=basePartition;
    code.row=row;
    code.col=col;
    isea->decode(code,coord);
    originLon=coord.lon;
    originLat=coord.lat;

    int gridTiffIndex=getGridTiffIndex(basePartition, row, col, griTiffCodes, griTiffCount, originLevel);
    GridTransform gridTrans=getGridTrans(griTiffCodes[gridTiffIndex], originLevel);
    T elev = getValueFromGridTIFF(row,col, eleValues+gridTiffIndex*GridTIFFSize(gridSize), gridTrans, gridSize);

    if(elev<minAltitude)
        elev=0;
    //格元空间坐标
    SpatialCoord shetCoord{0, 0, 0};
    LLE2XYZ(originLon, originLat, (double)elev, shetCoord.x, shetCoord.y, shetCoord.z);

    double sPitch = calcPitch(shetCoord.x, shetCoord.y, shetCoord.z, satSpaceCoord.x, satSpaceCoord.y, satSpaceCoord.z);
    if (sPitch < 0)
        return 0;

    //卫星经纬度坐标
    double sLon = 0, sLat = 0,sElev=0;
    XYZ2LLE(satSpaceCoord.x,satSpaceCoord.y,satSpaceCoord.z, sLon, sLat,sElev);

    //坐标方位角
    double yaw = getYaw(originLon, originLat, sLon, sLat);
    //总长
    double totalDistance = getDistance(originLon, originLat, sLon, sLat);

    totalDistance = totalDistance < (double)(MaxDistance) ? totalDistance : (double)(MaxDistance);
    //初始步长
    double cellStepDistance =1000 * earthRadius * 2 * M_PI / (5 * pow(2, level + 1));

    double sumStepDis = cellStepDistance;

    double vtan = tan(sPitch);


    double interElev=elev;
    while (sumStepDis <= totalDistance)
    {
        //计算间隔点经纬度信息
        double interLon = 0, interLat = 0;
        getNextCoord(originLon, originLat, yaw, sumStepDis, interLon, interLat);
        int interBasePartition=0,interRow=0,interCol=0;

        code.level=level;
        code.elementType=DGGSElementType::Cell;

        coord.lon=interLon;
        coord.lat=interLat;
        isea->encode(coord,code);
        interBasePartition=code.basePartition;
        interRow=code.row;
        interCol=code.col;

        //由间隔点经纬度获取高程
        gridTiffIndex=getGridTiffIndex(interBasePartition, interRow, interCol, griTiffCodes, griTiffCount, originLevel);

        if(gridTiffIndex!=-1)
        {
            gridTrans=getGridTrans(griTiffCodes[gridTiffIndex], originLevel);
            interElev= getValueFromGridTIFF(interRow, interCol, eleValues + gridTiffIndex * GridTIFFSize(gridSize), gridTrans, gridSize);
        }
       if(interElev<minAltitude)
            interElev=0.0;
        //间隔点空间坐标
        LLE2XYZ(interLon, interLat, interElev, shetCoord.x, shetCoord.y, shetCoord.z);

        //计算遮挡点与初始点的高度角
        double stepPitch = calcPitch(satSpaceCoord.x, satSpaceCoord.y, satSpaceCoord.z, shetCoord.x, shetCoord.y, shetCoord.z);
        if (stepPitch > sPitch)
            return 0;

        //下一次步长
        double dHeight = elev + sumStepDis * vtan - interElev + (sumStepDis * sumStepDis) / (2 * earthRadius*1000);
        double stepDistance = (dHeight / vtan) * 0.8;
        if (stepDistance < cellStepDistance)
            stepDistance = cellStepDistance;
        else if (stepDistance > MaxDistance/4)
            stepDistance = MaxDistance/4;
        sumStepDis += stepDistance;
    }
    return 1;
}

__device__ float calPdopValue(SpatialCoord* satCoords, int* satVis, const int satVisNum, SpatialCoord tagSpaceCoord)
{
    float pdopValue = 99;
    if (satVisNum >= 4)
    {
        double Q[satMaxCount * 4];
        for (int i = 0; i < satVisNum; i++)
        {
            double xSat = satCoords[satVis[i]].x;
            double ySat = satCoords[satVis[i]].y;
            double zSat = satCoords[satVis[i]].z;
            double u = xSat - tagSpaceCoord.x;
            double v = ySat - tagSpaceCoord.y;
            double w = zSat - tagSpaceCoord.z;
            Q[i * 4 + 0] = u / sqrt(u * u + v * v + w * w);
            Q[i * 4 + 1] = v / sqrt(u * u + v * v + w * w);
            Q[i * 4 + 2] = w / sqrt(u * u + v * v + w * w);
            Q[i * 4 + 3] = 1;
        }
        double QTQ[4*4];
        for (int i = 0; i < 16; ++i)
        {
            int ir = i / 4;
            int ic = i - ir * 4;
            double* pr = QTQ + ir * 4 + ic;
            pr[0] = 0;

            for (int j = 0; j < satVisNum; ++j)
                pr[0] += Q[j * 4 + ir] * Q[j * 4 + ic];
        }
        double H[4*4];
        int iret = getMatrixInverse(QTQ, H);

        if (iret == 1)
        {
            double trace = H[0] + H[5] + H[10];
            if (trace > 0.0)
                pdopValue = sqrt(float(trace));
        }
    }
    return pdopValue;

}

template<class T>
__global__ void getSatVisMat(SpatialCoord*satSpatialCoords,
                             int satIndex,
                             T* eleValues,
                             DGGSGridTIFFCode* griTIffCodes,
                             int griTiffCount,
                             int originLevel,
                             DGGSGridTIFFCode tagGridTiffCode,
                             int* visValues,
                             ICuDGGSKernel *isea
                             )
{
    int gridSize=(int)pow(2,originLevel);
    GridTransform tagTrans= getGridTrans(tagGridTiffCode,originLevel);
    //根据卫星的ID获取坐标
    SpatialCoord currentSat=satSpatialCoords[satIndex];

    for ( unsigned int i = blockDim.x * blockIdx.x + threadIdx.x; i < gridSize*gridSize; i += blockDim.x * gridDim.x)
    {
        unsigned int row = (unsigned int) (i / gridSize);
        unsigned int col = (unsigned int) (i - row * gridSize);

        int isVis=isVisible(currentSat, tagTrans.basePartition,row+tagTrans.rowOrigin,col+tagTrans.colOrigin, tagTrans.level, eleValues, griTIffCodes, griTiffCount, originLevel, isea);
        int x, y;
        code2XY(row, col, 0, 0, gridSize, x, y);
        if(x<0||x>gridSize*2-1||y<0||y>gridSize*2-1)
            printf("x=%d,y=%d is wrong\n",x,y);
        else {
            visValues[satIndex * gridSize * 2 * gridSize * 2 + (gridSize * 2 - 2 - y) * (gridSize * 2) + x] = isVis;
            visValues[satIndex * gridSize * 2 * gridSize * 2 + (gridSize * 2 - 2 - y) * (gridSize * 2) + x + 1] = isVis;
        }
    }
}

template<class T>
__global__ void getPDOPMat(T*eleValues,
                           int* visValues,
                           SpatialCoord*satSpatialCoords,
                           int satCount,
                           float*tagPDOPValues,
                           GridTransform tagTrans,
                           int gridSize,
                           ICuDGGSKernel* isea
                           )
{
    int basePartition=tagTrans.basePartition;
    for (unsigned int index = blockDim.x * blockIdx.x + threadIdx.x; index < gridSize * gridSize; index += blockDim.x * gridDim.x)
    {
        unsigned int row = (int) (index / gridSize);
        unsigned int col = (int) (index - row * gridSize);
        int level = tagTrans.level;

        //用于记录可见卫星
        int satVisNum = 0;
        int satVis[satMaxCount] = {0};

        int x, y;
        code2XY(row, col, 0, 0, gridSize, x, y);

        size_t  gridIndex=(gridSize * 2 - 2 - y) * (gridSize * 2) + x;

        for (int i = 0; i < satCount; i++) {
            int visible = visValues[i * gridSize * 2 * gridSize * 2 + gridIndex];
            if (visible == 1) {
                satVis[satVisNum] = i;
                satVisNum++;
            }
        }
        double lon=0, lat=0, elev=0;
        DGGSElementCode code;
        code.basePartition = basePartition;
        code.row = row;
        code.col = col;
        code.level = level;
        code.elementType=DGGSElementType::Cell;
        GeoCoord coord;
        isea->decode(code,coord);

        lon=coord.lon;
        lat=coord.lat;
        elev =(double) eleValues[gridIndex];

        if(elev<minAltitude)
            elev=0.0;

        SpatialCoord tagCoord{};
        LLE2XYZ(lon, lat, elev, tagCoord.x, tagCoord.y, tagCoord.z);

        float pdopValue = calPdopValue(satSpatialCoords, satVis, satVisNum, tagCoord);

        tagPDOPValues[gridIndex] = pdopValue;
        tagPDOPValues[gridIndex+1] = pdopValue;
    }
}

#pragma endregion

#pragma region buffer function Declaration

/// <summary>
/// 获取GridTIFF的空间边界
/// <summary>
///  <param name="gridTiffCode">GridTIFFCode</param>
/// <returns>DGGridDiamond</returns>
DGGridDiamond getSpaceBoundFromGridTiffCode(DGGSGridTIFFCode gridTiffCode);

/// <summary>
/// 获取缓冲区GridTIFF编码
/// </summary>
bool getGridTiffFromBound();

/// <summary>
/// 读取GridTIFF文件，加载缓冲区数据
/// </summary>
/// <param name="srcFileName">源文件</param>
/// <returns></returns>
bool loadFromFile(const char* srcFileName);

/// <summary>
/// 获取当前GridTIFF在缓冲区中的索引
/// </summary>
///  <returns>索引</returns>
int getTagGridTiffIndex();

#pragma endregion

#pragma region buffer function Implementation

DGGridDiamond getSpaceBoundFromGridTiffCode(DGGSGridTIFFCode gridTiffCode)
{
    DGGridDiamond diamond;
    DGGSElementCode code;
    code.elementType=pElementType;

    gridTiffCode2GridCode(gridTiffCode, code, globalOriginLevel);
    diamond.basePartition=gridTiffCode.basePartition;

    auto rowMin = code.row;
    auto colMin = code.col;
    long long int rowMax;
    long long int colMax;
    if (code.level < globalOriginLevel)
    {
        rowMax = (long long)std::pow(2,code.level);
        colMax = (long long) std::pow(2,code.level);
    }
    else
    {
        rowMax = rowMin + (long long)std::pow(2, globalOriginLevel);
        colMax = colMin + (long long)std::pow(2, globalOriginLevel);
    }


    auto codes=new DGGSElementCode[4];
    for(int i=0;i<4;i++)
    {
        codes[i].basePartition=gridTiffCode.basePartition;
        codes[i].level=code.level;
        codes[i].elementType=DGGSElementType::Cell;
    }

    codes[0].row=rowMin;
    codes[0].col= colMax - 1;

    codes[1].row= rowMax - 1;
    codes[1].col= colMin;

    codes[2].row= rowMin;
    codes[2].col= colMin;

    codes[3].row= rowMax - 1;
    codes[3].col= colMax - 1;

    auto pKernel= CreateKernelObj(pGridType);

    auto coord= new GeoCoord[4];
    for(int i=0;i<4;i++)
        pKernel->decode(codes[i],coord[i]);

    //上
    diamond.latMax = coord[0].lat;

    //下
    diamond.latMin = coord[1].lat;

    //左
    diamond.lonMin = coord[2].lon;

    //右
    diamond.lonMax = coord[3].lon;

    delete[]coord;
    delete[]codes;

    return diamond;
}

bool getGridTiffFromBound()
{
    //存储GridTIFF编码
    std::vector<DGGSGridTIFFCode> resultVector;

    double differ= calculateDiffer(pLevel);

    int xSize=0;
    if(bufferBound[3]-bufferBound[2]>0)
        xSize=std::ceil((bufferBound[3]-bufferBound[2])/differ);
    else
        xSize=std::ceil((bufferBound[3]+360-bufferBound[2])/differ);
    int ySize=std::ceil((bufferBound[0]-bufferBound[1])/differ);

    auto pKernel= CreateKernelObj(pGridType);

    int xStep= 360.0 / (pow(2, pLevel - globalOriginLevel)) / differ / 10;
    int yStep= 180.0 / (pow(2, pLevel - globalOriginLevel)) / differ / 10;

    for(int x=0;x<xSize+xStep;x+=xStep)
    {
        for (int y = 0; y < ySize+yStep; y+=yStep)
        {
            double lon = bufferBound[2] + x * differ;
            double lat = bufferBound[0] - y * differ;

            if(lon<-180)lon=360.0+lon;
            if(lon>180)lon=lon-360.0;

            DGGSElementCode tempCode;
            tempCode.level = pLevel;
            tempCode.elementType=pElementType;
            pKernel->encode(GeoCoord{lon,lat},tempCode);

            if(tempCode.basePartition==0||tempCode.basePartition==11)
                continue;

            DGGSGridTIFFCode tileCode;
            gridCode2GridTiffCode(tempCode, tileCode, globalOriginLevel);

            int i=0;
            for(;i<resultVector.size();i++)
                if(tileCode==resultVector[i])
                    break;

            if(i==resultVector.size())
                resultVector.push_back(tileCode);
        }
    }

    DestroyKernelObj(pKernel);

    //拷贝结果
    bufferGridTIFFCount = (int)resultVector.size();
    bufferGridTIFFCodes = new DGGSGridTIFFCode[bufferGridTIFFCount];
    for (int i = 0; i < bufferGridTIFFCount; i++)
    {
        bufferGridTIFFCodes[i] = resultVector[i];
    }

    return true;
}

bool loadFromFile(const char* srcFileName)
{
    //获取文件夹路径
    std::string srcFilePath =  fs::path(srcFileName).parent_path().string();

    //当前GridTIFF
    std::string gridCodeStr = fs::path(srcFileName).stem().string();
    globalGridCode=::getGridTiffCode(gridCodeStr);
    pLevel=globalGridCode.level;

    globalOriginLevel= globalGridCode.level - (gridCodeStr.length() - 2);
    globalGridSize =pow(2, globalOriginLevel);
    //获取缓冲区范围
    DGGridDiamond diamond=getSpaceBoundFromGridTiffCode(globalGridCode);
    diamond.latMax+=BufferSize*radToDeg;
    diamond.lonMin-=BufferSize*radToDeg;
    diamond.latMax+=BufferSize*radToDeg;
    diamond.latMin-=BufferSize*radToDeg;
    if(diamond.latMax>90)diamond.latMax=90.0;
    if(diamond.latMin<-90)diamond.latMin=-90.0;

    //上
    bufferBound[0]=diamond.latMax;
    //下
    bufferBound[1]=diamond.latMin;
    //左
    bufferBound[2]=diamond.lonMin;
    //右
    bufferBound[3]=diamond.lonMax;
    getGridTiffFromBound();

    size_t left = std::string(srcFileName).find_last_of(".");
    std::string exFormat = std::string(srcFileName).substr(left);

    bool  isMalloc =false;

    for(int i=0; i < bufferGridTIFFCount; i++)
    {
        auto currentTiffCode=bufferGridTIFFCodes[i];
        std::string currentFileName= srcFilePath + "/" + getGridTiffString(currentTiffCode, globalOriginLevel) + exFormat;

        if(!fs::exists(currentFileName.c_str()))
            continue;

        GridTIFF gridTiff;
        gridTiff.loadFromFile(currentFileName.c_str());

        //为存储数据分配缓冲区
        if(!isMalloc)
        {
            globalDataType = gridTiff.getDataType();
            size_t allDataSize = GDALGetDataTypeSizeBytes(static_cast<GDALDataType>(globalDataType)) * GridTIFFSize(globalGridSize) * bufferGridTIFFCount;
            bufferGridTIFFValues= CPLMalloc(allDataSize);
            if (bufferGridTIFFValues != nullptr) {
                memset(bufferGridTIFFValues, 0, allDataSize);
            }
            isMalloc=true;
        }

        size_t dataSize = GDALGetDataTypeSizeBytes(static_cast<GDALDataType>(globalDataType)) * GridTIFFSize(globalGridSize);
        void*buffer=gridTiff.getBandValue(0);
        memcpy((char*)bufferGridTIFFValues + dataSize * i, buffer, dataSize);
    }
    return true;
}

int getTagGridTiffIndex()
{
    for(int i=0; i < bufferGridTIFFCount; i++)
        if(bufferGridTIFFCodes[i] == globalGridCode)
            return i;
    return -1;
}

#pragma endregion

/// <summary>
///计算单个GridTIFF的PDOP
/// </summary>
/// <param name="srcFileName">源文件</param>
/// <param name="tagFileName">目标文件</param>
/// <pram name="satSpaceCoords">卫星空间坐标</param>
/// <pram name="satCount">卫星数量</param>
/// <returns></returns>
template<class T>
void calSingleGridTiffPdop(const char* srcFileName, const char*tagFileName, SpatialCoord *satSpaceCoords, int satCount)
{
    loadFromFile(srcFileName);

    //GridTIFF数据
    int gridSize=globalGridSize;
    int bufferGridCount=bufferGridTIFFCount;
    T* bufferValues=reinterpret_cast<T*>(bufferGridTIFFValues);
    T* bufferValues_d=nullptr;
    CHECK(cudaMalloc((T**)&bufferValues_d, sizeof(T) * bufferGridCount*GridTIFFSize(gridSize)));
    CHECK(cudaMemcpy((T**)bufferValues_d, bufferValues, GridTIFFSize(gridSize) * sizeof(T), cudaMemcpyHostToDevice));

    //可见性数据
    int* visValues_d=nullptr;
    CHECK(cudaMalloc((int**)&visValues_d, sizeof(int) * satCount * GridTIFFSize(gridSize) ));
    cudaMemset((int**)visValues_d,0,GridTIFFSize(gridSize) * satCount * sizeof(int));

    //相关GridTIFFCode
    auto gridTiffCodes=bufferGridTIFFCodes;
    DGGSGridTIFFCode tagGridTiffCode=globalGridCode;
    int originLevel=globalOriginLevel;
    DGGSGridTIFFCode *gridTiffCodes_d= nullptr;
    CHECK(cudaMalloc((DGGSGridTIFFCode**)&gridTiffCodes_d, sizeof(DGGSGridTIFFCode) * bufferGridCount));
    CHECK(cudaMemcpy(gridTiffCodes_d, gridTiffCodes, sizeof(DGGSGridTIFFCode) * bufferGridCount, cudaMemcpyHostToDevice));

    //卫星坐标
    SpatialCoord *satSpaceCoords_d= nullptr;
    CHECK(cudaMalloc((SpatialCoord**)&satSpaceCoords_d, sizeof(SpatialCoord) * satCount));
    CHECK(cudaMemcpy(satSpaceCoords_d, satSpaceCoords, sizeof(SpatialCoord) * satCount, cudaMemcpyHostToDevice));

    //PDOP值
    int tagIndex=getTagGridTiffIndex();
    float* tagPDOPValues_d=nullptr;
    CHECK(cudaMalloc((float**)&tagPDOPValues_d, sizeof(float) * GridTIFFSize(gridSize)));

    //获取GridTransform
    GridTransform tagTrans{};
    DGGSElementCode tempCode;
    gridTiffCode2GridCode(tagGridTiffCode, tempCode, globalOriginLevel);
    tagTrans.basePartition=tagGridTiffCode.basePartition;
    tagTrans.level=tagGridTiffCode.level;
    tagTrans.rowOrigin=(int)tempCode.row;
    tagTrans.colOrigin=(int)tempCode.col;

    T* eleValue=bufferValues+tagIndex*GridTIFFSize(gridSize);
    T*eleValue_d=nullptr;
    CHECK(cudaMalloc((T**)&eleValue_d, sizeof(T) * GridTIFFSize(gridSize)));

    //计算可见性
    for(int satIndex=0;satIndex<satCount;satIndex++)
    {
        getSatVisMat<<<pGridDimX, pBlockDimX>>>(
                satSpaceCoords_d, satIndex, bufferValues_d, gridTiffCodes_d,bufferGridCount, originLevel, tagGridTiffCode, visValues_d, pCuDGGSKernel);
        cudaDeviceSynchronize();
    }

    //设置GridTIFF
    GridTIFF gridTiff;
    if (!(std::filesystem::exists(tagFileName)))
    {
        gridTiff.setSize(gridSize);
        gridTiff.setDataType(DGGSDataType::DGGS_Float32);
        gridTiff.setBandCount(1);
        gridTiff.setTrans(tagTrans);
        gridTiff.setNoValue(noData);
        gridTiff.setElementType(pElementType);
        gridTiff.setGridType(pGridType);
    }
    else
    {
        gridTiff.loadFromFile(tagFileName);
    }
    auto pdopValues=reinterpret_cast<T*>(gridTiff.getBandValue(0));
    std::fill(pdopValues, pdopValues + GridTIFFSize(globalGridSize), noData);

    CHECK(cudaMemcpy(tagPDOPValues_d, pdopValues, sizeof(float) * GridTIFFSize(gridSize), cudaMemcpyHostToDevice));
    CHECK(cudaMemcpy(eleValue_d, eleValue, sizeof(T) * GridTIFFSize(gridSize), cudaMemcpyHostToDevice));

    //计算PDOP
    getPDOPMat<<<pGridDimX, pBlockDimX>>>(eleValue_d,visValues_d, satSpaceCoords_d,satCount, tagPDOPValues_d,tagTrans,gridSize,pCuDGGSKernel);
    cudaDeviceSynchronize();

    CHECK(cudaMemcpy(pdopValues, tagPDOPValues_d, sizeof(float) * GridTIFFSize(gridSize), cudaMemcpyDeviceToHost));

    gridTiff.saveToFile(tagFileName);

    //释放内存
    cudaFree(bufferValues_d);
    cudaFree(visValues_d);
    cudaFree(gridTiffCodes_d);
    cudaFree(satSpaceCoords_d);
    cudaFree(tagPDOPValues_d);
    cudaFree(eleValue_d);

    CPLFree(bufferGridTIFFValues);
    delete[]bufferGridTIFFCodes;
}


/// PDOP计算
void calculatePdop(const char* srcGridFilePath,
                   const char* gridExtFormat,
                   const char* srcSatFilePath,
                   const char* satExtFormat,
                   const char* satSyetem,
                   double timeStep,
                   const char* tagFilePath,
                   int gridDimX,
                   int blockDimX
                   )
{

    //解析星历数据并存入temp
    fs::path exeDir = fs::current_path();
    std::string tempDir = exeDir.string()+"/temp";
    fs::create_directories(tempDir);
    std::string* sp3FileNames;
    int sp3fileCount;
    getAllFilesFromPath(srcSatFilePath,satExtFormat,sp3FileNames,sp3fileCount);
    for(int i=0;i<sp3fileCount;i++)
        parseSp3File(sp3FileNames[i].c_str(),tempDir.c_str(),satSyetem,timeStep);
    delete[]sp3FileNames;

    //获取所有GridTiff文件
    std::string * gridTiffFiles= nullptr;
    int gridFileCount=0;
    getAllFilesFromPath(srcGridFilePath, gridExtFormat, gridTiffFiles, gridFileCount);

    //获取卫星数据
    std::string * satFiles= nullptr;
    int satFileCount=0;
    getAllFilesFromPath(tempDir, "", satFiles, satFileCount);

    //获取数据类型、格网类型等参数
    DGGSDataType dataType;
    {
        GridTIFF gridTiff;
        gridTiff.loadFromFile(gridTiffFiles[0].c_str());
        dataType=gridTiff.getDataType();
        pGridType = gridTiff.getGridType();
        pElementType =gridTiff.getElementType();
        pCuDGGSKernel= CreateCuKernelObj(pGridType);
    }

    //输入参数
    pGridDimX = gridDimX;
    pBlockDimX = blockDimX;


    //遍历所有星历数据
    for(int satFileIndex=0;satFileIndex<satFileCount;satFileIndex++)
    {
        std::cout<<satFileIndex<<" "<<satFiles[satFileIndex]<<std::endl;
        //星历文件名
        size_t satLeft = satFiles[satFileIndex].find_last_of("/\\");
        std::string satFileName = (satLeft != std::string::npos)?satFiles[satFileIndex].substr(satLeft + 1): satFiles[satFileIndex];

        //输出文件路径
        std::string currentTagFilePath=std::string(tagFilePath)+"/"+satFileName;
        if (!std::filesystem::exists(currentTagFilePath))
            std::filesystem::create_directories(currentTagFilePath);

        //加载星历文件
        SpatialCoord *satSpaceCoords= nullptr;
        int satCount=0;
        loadSatFile(satFiles[satFileIndex].c_str(),satSpaceCoords,satCount);

#define HANDLE_DATA_TYPE(T, D) \
        case  D:\
        {\
            calSingleGridTiffPdop<T>(srcFileName.c_str(),tagFileName.c_str(),satSpaceCoords,satCount);                   \
            break;                         \
        }

        //逐个GridTIFF计算PDOP
        for(int tiffIndex=0;tiffIndex<gridFileCount;tiffIndex++)
        {
            size_t left = gridTiffFiles[tiffIndex].find_last_of("/\\");
            size_t right = gridTiffFiles[tiffIndex].find_last_of('.');
            std::string fileName=(right != std::string::npos) ? gridTiffFiles[tiffIndex].substr(left+1,right-left-1) : gridTiffFiles[tiffIndex];
            std::string srcFileName=std::string(srcGridFilePath)+"/"+fileName+gridExtFormat;
            std::string tagFileName=std::string(tagFilePath)+"/"+satFileName+"/"+fileName+gridExtFormat;
            switch (dataType) {
                HANDLE_DATA_TYPE(unsigned short, DGGS_Byte)
                HANDLE_DATA_TYPE(short, DGGS_Int16)
                HANDLE_DATA_TYPE(unsigned int, DGGS_UInt32)
                HANDLE_DATA_TYPE(int, DGGS_Int32)
                HANDLE_DATA_TYPE(float, DGGS_Float32)
                HANDLE_DATA_TYPE(double, DGGS_Float64)
                default:
                    break;
            }
        }
    }
    delete[]satFiles;
    delete[]gridTiffFiles;
    DestroyCuKernelObj(pCuDGGSKernel);
}